19  CJS [t,t] - growth model

Simple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to compare with Neural Network CMR/growth model

In all the models below, 1 = not observed and 2 = observed in the input encounter data.
Encounter data are available here in the eh.csv file. Weight data are in weight.csv

19.0.1 modelNum 1: g(i,t) model

Growth rate (gr) model structure:

gr(t,i) <- grInt(t,i)

Model code is in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_makeFile.R
Not using this: Targets are loaded in R/modelCMR_tt_growth_NN_OB.R and saved as modelCMR_tt_growth_NN_OB_target

Model runs:

19.0.1.1 Retrieve model results

Code
library(targets)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
modelNum = 1 # growth only
dummy=2

#load('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_mostRecent.RData')
load(paste0('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_model', modelNum, '_mostRecent.RData'))
out_NN_OBmod1 <- d
rm(d)

#Input data
eh <- tar_read(eh_OB_2002_2014_target)

19.0.1.2 Model code

In the model code, a value of 1 for z or in gamma or omega signifies the individual is alive and a value of 2 signifies the individual is dead. y[i,j] is the observed data (encounter history file).

Code
out_NN_OBmod1$modelCode
{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AIS[first[i]], 
            sd = 2)
        weightDATA[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = 0.2)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATA[i, t + 1] ~ dnorm(weight[i, t + 1], sd = 0.2)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ dnorm(0, sd = 1000)
    }
}

19.0.1.3 Model statistics

Run statistics
nIter nBurnin nChains thinRate
4000 2000 2 5

[1] “Run time = 32.138 mins”

19.0.1.4 Plot traces

Code
  priors <- rnorm(out_NN_OBmod1$runData$nIter * out_NN_OBmod1$runData$nChains, 0, 1000)
  MCMCtrace(object = out_NN_OBmod1$mcmc,
            #ISB = FALSE,
            #exact = TRUE, 
            params = c("grIntT"),
            pdf = FALSE, 
            priors = priors)

19.0.1.5 Summary data

Code
  MCMCplot(object = out_NN_OBmod1$mcmc, params = c("grIntT"))

Summary values

Code
#|
  (summaryMod1_tt_growth <- MCMCsummary(object = out_NN_OBmod1$mcmc, params = c("grIntT"), round = 3) %>%
    rownames_to_column(var = "var"))
          var   mean    sd   2.5%    50%  97.5% Rhat n.eff
1   grIntT[1]  1.169 0.028  1.119  1.168  1.227 1.01   197
2   grIntT[2]  1.122 0.022  1.080  1.123  1.165 1.01   168
3   grIntT[3]  9.277 0.028  9.222  9.278  9.329 1.01    91
4   grIntT[4]  3.327 0.023  3.282  3.327  3.374 1.03   212
5   grIntT[5]  0.584 0.033  0.521  0.585  0.647 1.05   135
6   grIntT[6]  1.557 0.036  1.485  1.556  1.622 1.03   134
7   grIntT[7] 17.120 0.045 17.030 17.121 17.212 1.13    59
8   grIntT[8]  7.398 0.059  7.288  7.397  7.509 1.00    86
9   grIntT[9] -1.513 0.080 -1.665 -1.512 -1.333 1.05    39
10 grIntT[10]  4.099 0.083  3.942  4.096  4.257 1.13    43
11 grIntT[11] 22.872 0.120 22.625 22.865 23.101 1.10    21
Code
  # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
  source('./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R')
  
  #summaryMod1_tt_growth_z0 <- MCMCSummaryRMNA(object = out_NN_OBmod1$mcmc, params = c("weight"), round = 3) %>%
  summaryMod1_tt_growth_zGR <- MCMCSummaryRMNA(object = out_NN_OBmod1$mcmc, params = c("gr"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod1_tt_growth_zW <- MCMCSummaryRMNA(object = out_NN_OBmod1$mcmc, params = c("weight"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod1_tt_growth_z0 <- summaryMod1_tt_growth_zW |> 
    left_join(summaryMod1_tt_growth_zGR, by = c("ind", "t"), suffix = c("_weight", "_gr"))

Add other variables to summary values

Code
  ehLong <- 
    eh$eh |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      enc = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))

  firstLong <- eh$first |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(first = value)
  
  lastLong <- eh$last |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(last = value)
    
  cohortLong <- eh$cohort |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) 
  
  weightLong <- 
    eh$weight |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      weightDATA = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  summaryMod1_tt_growth_z <- summaryMod1_tt_growth_z0 |> 
    left_join(ehLong) |> 
    mutate(
      # meanM1 = mean - 1,
      # pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
      enc8 = ifelse(enc == 1, 0.8, 0)
    ) |> 
    left_join(firstLong) |> 
    left_join(lastLong) |> 
    left_join(cohortLong) |> 
    left_join(weightLong) |> 
    mutate(
      resid = weightDATA - mean_weight
    ) 
Code
  resids <- summaryMod1_tt_growth_z |> 
    group_by(ind) |> 
    summarize(
      meanResid = mean(abs(resid), na.rm = TRUE),
      n = n()
    ) |> 
    ungroup()

  summaryMod1_tt_growth_zR <- 
    summaryMod1_tt_growth_z |> 
    left_join(resids)
Code
#  write.csv(summaryMod1_tt_growth_zR, './data/outForObservable/summaryMod1_tt_growth_zR.csv')
Code
ojs_define(numTags_OJS = dim(eh$tags)[1])
ojs_define(summaryMod1_tt_growth_zR_OJS = (summaryMod1_tt_growth_zR))

19.0.1.6 Plot predicted and observed masses for selected individuals

Select one or more individuals:

Code
viewof selectInd = Inputs.select([...Array(numTags_OJS).keys()], {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})

summaryMod1_tt_growth_zR_OJS_T = transpose(summaryMod1_tt_growth_zR_OJS)
summaryMod1_tt_growth_zR_OJS_T_selected = summaryMod1_tt_growth_zR_OJS_T.filter((d) =>
  selectInd.includes(d.ind)
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort. Green dots are estimated growth rates.

Code
Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Body mass (mg)" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    Plot.line(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    Plot.dot(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "t",
      y: "mean_gr",
      fill: "green"
    }),
    Plot.dot(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "t",
      y: "weightDATA",
      fill: "orange"
    }),
    Plot.ruleX(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "first",
      y: 3,
      stroke: "green"
    }),
    Plot.ruleX(summaryMod1_tt_growth_zR_OJS_T_selected, {
      x: "last",
      y: 3,
      stroke: "red"
    }),
    Plot.text(summaryMod1_tt_growth_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 60,
        frameAnchor: "top-left",
        text: (d) => "Cohort = " + d.cohort
      })
    ),
    Plot.text(summaryMod1_tt_growth_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 56,
        frameAnchor: "top-left",
        text: (d) => "Residual = " + d.meanResid
      })
    )
  ],
  facet: {
    data: summaryMod1_tt_growth_zR_OJS_T_selected,
    x: "ind"
  }
})

19.0.2 modelNum 2: phi(i,t) + g(i,t), p(i,t) model

Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)

Probability of survival (phi) model structure:

phi(t,i) <- phiInt(t,i)

Probability of capture (p) model structure:

p(t,i) <- pInt(t,i)

Growth rate (gr) model structure:

gr(t,i) <- grInt(t,i)

Survival/growth rate interaction:

Additive

Model code is in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_makeFile.R
Not using this: Targets are loaded in R/modelCMR_tt_growth_NN_OB.R and saved as modelCMR_tt_growth_NN_OB_target

Model runs:

19.0.2.1 Retrieve model results

Code
library(targets)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
modelNum = 2 # phi + growth

#load('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_mostRecent.RData')
load(paste0('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_model', modelNum, '_mostRecent.RData'))
out_NN_OBmod2 <- d
rm(d)

#Input data
eh <- tar_read(eh_OB_2002_2014_target)

19.0.2.2 Model code

In the model code, a value of 1 for z or in gamma or omega signifies the individual is alive and a value of 2 signifies the individual is dead. y[i,j] is the observed data (encounter history file).

Code
out_NN_OBmod2$modelCode
{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(3.7, sd = 1.7)
        weightDATA[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = 0.2)
        for (t in first[i]:(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]
            weightDATA[i, t + 1] ~ dnorm(weight[i, t + 1], sd = 0.2)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] ~ dnorm(grIntT[t], sd = 0.33)
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ T(dnorm(0, sd = 1.5), -3.5, 3.5)
        grIntTOut[t] <- ilogit(grIntT[t])
    }
    delta[1] <- 1
    delta[2] <- 0
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gamma[1, 1, t, i] <- phi[i, t]
            gamma[1, 2, t, i] <- 1 - phi[i, t]
            gamma[2, 1, t, i] <- 0
            gamma[2, 2, t, i] <- 1
        }
    }
    for (t in 1:(T - 1)) {
        p[t] ~ dunif(0, 1)
        omega[1, 1, t] <- 1 - p[t]
        omega[1, 2, t] <- p[t]
        omega[2, 1, t] <- 1
        omega[2, 2, t] <- 0
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            logit(phi[i, t]) <- phiInt[i, t]
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
        }
    }
    for (t in 1:(T - 1)) {
        phiIntT[t] ~ T(dnorm(0, sd = 1.5), -3.5, 3.5)
        phiIntTOut[t] <- ilogit(phiIntT[t])
    }
    for (i in 1:N) {
        z[i, first[i]] ~ dcat(delta[1:2])
        for (j in first[i]:(last[i] - 1)) {
            z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
            yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
        }
    }
}

19.0.2.3 Model statistics

Run statistics
nIter nBurnin nChains thinRate
5000 2000 2 5

[1] “Run time = 3.472 hours”

19.0.2.4 Plot traces

Code
  priors <- rnorm(out_NN_OBmod2$runData$nIter * out_NN_OBmod2$runData$nChains, 0, 1000)
  MCMCtrace(object = out_NN_OBmod2$mcmc,
            #ISB = FALSE,
            #exact = TRUE, 
            params = c("phiIntTOut", "grIntT"),
            pdf = FALSE, 
            priors = priors)

19.0.2.5 Summary data

Code
  MCMCplot(object = out_NN_OBmod2$mcmc, params = c("phiIntTOut"))

Code
  MCMCplot(object = out_NN_OBmod2$mcmc, params = c("p"))

Code
  MCMCplot(object = out_NN_OBmod2$mcmc, params = c("grIntT"))

Summary values

Code
#|
(summaryMod2_tt_growth <- MCMCsummary(object = out_NN_OBmod2$mcmc, params = c("phiIntTOut", "p", "grIntT"), round = 3) %>%
    rownames_to_column(var = "var"))
              var   mean    sd   2.5%    50%  97.5% Rhat n.eff
1   phiIntTOut[1]  0.914 0.015  0.883  0.913  0.943 1.11    60
2   phiIntTOut[2]  0.906 0.020  0.867  0.906  0.944 1.59    20
3   phiIntTOut[3]  0.855 0.020  0.814  0.855  0.897 1.26    38
4   phiIntTOut[4]  0.727 0.027  0.667  0.730  0.774 1.07    54
5   phiIntTOut[5]  0.714 0.036  0.645  0.715  0.787 1.10    26
6   phiIntTOut[6]  0.812 0.039  0.732  0.810  0.880 1.41    15
7   phiIntTOut[7]  0.904 0.037  0.827  0.918  0.954 1.36     7
8   phiIntTOut[8]  0.262 0.048  0.179  0.257  0.360 1.23    18
9   phiIntTOut[9]  0.680 0.085  0.497  0.687  0.834 1.03    11
10 phiIntTOut[10]  0.705 0.063  0.584  0.697  0.820 1.82    16
11 phiIntTOut[11]  0.670 0.255  0.148  0.741  0.948 8.00     3
12           p[1]  0.596 0.021  0.557  0.596  0.638 1.04   210
13           p[2]  0.448 0.019  0.413  0.448  0.486 1.16   171
14           p[3]  0.707 0.018  0.669  0.709  0.739 1.01   283
15           p[4]  0.658 0.021  0.615  0.658  0.699 1.01   200
16           p[5]  0.639 0.027  0.588  0.639  0.690 1.04   142
17           p[6]  0.492 0.029  0.438  0.492  0.552 1.02   124
18           p[7]  0.707 0.038  0.635  0.706  0.778 1.00    26
19           p[8]  0.770 0.052  0.666  0.774  0.866 1.04   128
20           p[9]  0.634 0.067  0.499  0.637  0.757 1.02   382
21          p[10]  0.660 0.083  0.502  0.660  0.816 1.21   151
22          p[11]  0.743 0.160  0.459  0.729  0.992 5.26    56
23      grIntT[1]  1.170 0.027  1.117  1.169  1.221 1.00   177
24      grIntT[2]  1.122 0.022  1.080  1.122  1.166 1.02   222
25      grIntT[3]  9.274 0.027  9.223  9.273  9.328 1.01    89
26      grIntT[4]  3.328 0.024  3.282  3.328  3.377 1.09   174
27      grIntT[5]  0.578 0.033  0.514  0.578  0.643 1.07    90
28      grIntT[6]  1.563 0.032  1.503  1.562  1.623 1.02   196
29      grIntT[7] 17.102 0.047 17.007 17.103 17.188 1.26    77
30      grIntT[8]  7.420 0.056  7.311  7.419  7.528 1.07    75
31      grIntT[9] -1.546 0.080 -1.695 -1.548 -1.399 1.21    55
32     grIntT[10]  4.091 0.089  3.908  4.088  4.284 1.11    41
33     grIntT[11] 22.853 0.151 22.590 22.836 23.185 1.74    15
Code
  # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
  source('./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R')
  
  #summaryMod2_tt_growth_z0 <- MCMCSummaryRMNA(object = out_NN_OBmod2$mcmc, params = c("weight"), round = 3) %>%
  summaryMod2_tt_w <- MCMCSummaryRMNA(object = out_NN_OBmod2$mcmc, params = c("weight"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod2_tt_z <- MCMCSummaryRMNA(object = out_NN_OBmod2$mcmc, params = c("z"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod2_tt_z0 <- summaryMod2_tt_z |> 
    left_join(summaryMod2_tt_w, by = c("ind", "t"), suffix = c("_z", "_weight"))

Add other variables to summary values

Code
  ehLong <- 
    eh$eh |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      enc = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  firstLong <- eh$first |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(first = value)
  
  lastLong <- eh$last |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(last = value)
    
  cohortLong <- eh$cohort |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) 
  
  weightLong <- 
    eh$weight |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      weightDATA = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  summaryMod2_tt_z <- summaryMod2_tt_z0 |> 
    left_join(ehLong) |> 
    mutate(
      # meanM1 = mean - 1,
      # pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
      enc8 = ifelse(enc == 1, 0.8, 0)
    ) |> 
    left_join(firstLong) |> 
    left_join(lastLong) |> 
    left_join(cohortLong) |> 
    left_join(weightLong) |> 
    mutate(
      resid = weightDATA - mean_weight
    ) 
    
  # summaryMod2_tt_z <- summaryMod2_tt_z0 |> 
  #   left_join(ehLong) |> 
  #   mutate(
  #     # meanM1 = mean - 1,
  #     # pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
  #     enc8 = ifelse(enc == 1, 0.8, 0)
  #   ) |> 
  #   left_join(firstLong) |> 
  #   left_join(lastLong) |> 
  #   left_join(cohortLong) 
Code
  resids_mod2 <- summaryMod2_tt_z |> 
    group_by(ind) |> 
    summarize(
      meanResid = mean(abs(resid), na.rm = TRUE),
      n = n()
    ) |> 
    ungroup()

  summaryMod2_tt_zR <- 
    summaryMod2_tt_z |> 
    left_join(resids_mod2)
Code
ojs_define(numTags_OJS_mod2 = dim(eh$tags)[1])
ojs_define(summaryMod2_tt_zR_OJS = (summaryMod2_tt_zR))
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod2))

19.0.2.6 Plot predicted and observed masses for selected individuals

Select one or more individuals:

Code
viewof selectInd_mod2 = Inputs.select([...Array(numTags_OJS).keys()], {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})

summaryMod2_tt_zR_OJS_T = transpose(summaryMod2_tt_zR_OJS)
summaryMod2_tt_zR_OJS_T_selected = summaryMod2_tt_zR_OJS_T.filter((d) =>
  selectInd_mod2.includes(d.ind)
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code
Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Body mass (mg)" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    Plot.line(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    <!-- Plot.dot(summaryMod2_tt_zR_OJS_T_selected, { -->
    <!--   x: "t", -->
    <!--   y: "mean_gr", -->
    <!--   fill: "green" -->
    <!-- }), -->
    Plot.dot(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "weightDATA",
      fill: "orange"
    }),
    Plot.ruleX(summaryMod2_tt_zR_OJS_T_selected, {
      x: "first",
      y: 3,
      "stroke": "green"
    }),
    Plot.ruleX(summaryMod2_tt_zR_OJS_T_selected, {
      x: "last",
      y: 3,
      "stroke": "red"
    }),
    Plot.text(summaryMod2_tt_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 60,
        frameAnchor: "top-left",
        text: (d) => "Cohort = " + d.cohort
      })
    ),
    Plot.text(summaryMod2_tt_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 56,
        frameAnchor: "top-left",
        text: (d) => "Residual = " + d.meanResid
      })
    )
  ],
  facet: {
    data: summaryMod2_tt_zR_OJS_T_selected,
    x: "ind"
  }
})

19.0.2.7 Plot survival

Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code
Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Probability of survival" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_z"
    }),
    Plot.dot(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "enc8",
      fill: "orange"
    }),
    Plot.line(summaryMod2_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_z"
    }),
    Plot.ruleX(summaryMod2_tt_zR_OJS_T_selected, {
      x: "first",
      y: 1,
      "stroke": "green"
    }),
    Plot.ruleX(summaryMod2_tt_zR_OJS_T_selected, {
      x: "last",
      y: 1,
      "stroke": "red"
    })
    ,
    Plot.text(summaryMod2_tt_zR_OJS_T_selected, 
      Plot.selectFirst({
        x: 11,
        y: 1,
        text: d => d.cohort
      })
    )
  ],
  facet: {
    data: summaryMod2_tt_zR_OJS_T_selected,
    x: "ind"
  }
})

19.0.3 modelNum 3: phi(i,t) * g(i,t), p(i,t) model

Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)

Probability of survival (phi) model structure:

phi(t,i) <- phiInt(t,i) + weight(t,i) + weight(t,i)^2

Probability of capture (p) model structure:

p(t,i) <- pInt(t,i)

Growth rate (gr) model structure:

gr(t,i) <- grInt(t,i)

Survival/growth rate interaction:

Multiplicative

Model code is in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_makeFile.R
Not using this: Targets are loaded in R/modelCMR_tt_growth_NN_OB.R and saved as modelCMR_tt_growth_NN_OB_target

Model runs:

19.0.3.1 Retrieve model results

Code
library(targets)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
modelNum = 3 # phi + growth

#load('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_mostRecent.RData')
load(paste0('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_model', modelNum, '_mostRecent.RData'))
out_NN_OBmod3 <- d
rm(d)

#Input data
eh <- tar_read(eh_OB_2002_2014_target)

19.0.3.2 Model code

In the model code, a value of 1 for z or in gamma or omega signifies the individual is alive and a value of 2 signifies the individual is dead. y[i,j] is the observed data (encounter history file).

Code
out_NN_OBmod3$modelCode
{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AIS[first[i]], 
            sd = 2)
        weightDATA[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = 0.2)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATA[i, t + 1] ~ dnorm(weight[i, t + 1], sd = 0.2)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ dnorm(0, sd = 1000)
    }
    delta[1] <- 1
    delta[2] <- 0
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gamma[1, 1, t, i] <- phi[i, t]
            gamma[1, 2, t, i] <- 1 - phi[i, t]
            gamma[2, 1, t, i] <- 0
            gamma[2, 2, t, i] <- 1
        }
    }
    for (t in 1:(T - 1)) {
        p[t] ~ dunif(0, 1)
        omega[1, 1, t] <- 1 - p[t]
        omega[1, 2, t] <- p[t]
        omega[2, 1, t] <- 1
        omega[2, 2, t] <- 0
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i, 
                t] * weightDATA[i, t] + phiBeta[2, i, t] * weightDATA[i, 
                t] * weightDATA[i, t]
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
            for (b in 1:2) {
                phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
            }
        }
    }
    for (t in 1:(T - 1)) {
        phiIntT[t] ~ T(dnorm(0, sd = 1.5), -3.5, 3.5)
        phiIntTOut[t] <- ilogit(phiIntT[t])
        phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
        phiBetaT[2, t] ~ dnorm(0, sd = 1/0.667)
    }
    for (i in 1:N) {
        z[i, first[i]] ~ dcat(delta[1:2])
        for (j in first[i]:(last[i] - 1)) {
            z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
            yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
        }
    }
}

19.0.3.3 Model statistics

Run statistics
nIter nBurnin nChains thinRate
5000 2500 2 5

[1] “Run time = 2.309 hours”

19.0.3.4 Plot traces

Code
  priors <- rnorm(out_NN_OBmod3$runData$nIter * out_NN_OBmod3$runData$nChains, 0, 1000)
  MCMCtrace(object = out_NN_OBmod3$mcmc,
            #ISB = FALSE,
            #exact = TRUE, 
            params = c("phiIntTOut", "grIntT"),
            pdf = FALSE, 
            priors = priors)

19.0.3.5 Summary data

Code
  MCMCplot(object = out_NN_OBmod3$mcmc, params = c("phiIntTOut"))

Code
  MCMCplot(object = out_NN_OBmod3$mcmc, params = c("phiBetaT"))

Code
  MCMCplot(object = out_NN_OBmod3$mcmc, params = c("p"))

Code
  MCMCplot(object = out_NN_OBmod3$mcmc, params = c("grIntT"))

Summary values

Code
#|
(summaryMod3_tt_growth <- MCMCsummary(object = out_NN_OBmod3$mcmc, params = c("phiIntTOut", "p", "grIntT"), round = 3) %>%
    rownames_to_column(var = "var"))
              var   mean    sd   2.5%    50%  97.5% Rhat n.eff
1   phiIntTOut[1]  0.830 0.043  0.739  0.831  0.905 1.96    27
2   phiIntTOut[2]  0.836 0.046  0.741  0.833  0.918 1.30    16
3   phiIntTOut[3]  0.825 0.037  0.738  0.827  0.892 1.21    19
4   phiIntTOut[4]  0.746 0.034  0.685  0.744  0.818 2.89    34
5   phiIntTOut[5]  0.757 0.038  0.684  0.760  0.822 1.01    29
6   phiIntTOut[6]  0.783 0.040  0.708  0.781  0.860 1.50    22
7   phiIntTOut[7]  0.862 0.031  0.783  0.870  0.907 1.56    28
8   phiIntTOut[8]  0.492 0.100  0.340  0.483  0.717 3.90    15
9   phiIntTOut[9]  0.544 0.108  0.344  0.528  0.752 1.21     7
10 phiIntTOut[10]  0.643 0.139  0.371  0.645  0.832 4.30    14
11 phiIntTOut[11]  0.576 0.122  0.373  0.560  0.821 1.45    10
12           p[1]  0.592 0.020  0.554  0.592  0.631 1.00   244
13           p[2]  0.408 0.018  0.374  0.408  0.444 1.00    98
14           p[3]  0.704 0.019  0.666  0.704  0.740 1.02   211
15           p[4]  0.656 0.021  0.615  0.655  0.696 1.15   162
16           p[5]  0.629 0.026  0.578  0.630  0.679 1.02    95
17           p[6]  0.504 0.026  0.454  0.505  0.555 1.10   192
18           p[7]  0.747 0.032  0.682  0.747  0.809 1.00   158
19           p[8]  0.783 0.046  0.684  0.785  0.866 1.02   289
20           p[9]  0.609 0.071  0.463  0.611  0.737 1.01   173
21          p[10]  0.669 0.086  0.501  0.668  0.833 1.00   120
22          p[11]  0.763 0.160  0.432  0.786  0.994 2.76    69
23      grIntT[1]  0.114 0.029  0.057  0.116  0.162 1.13   128
24      grIntT[2]  0.133 0.025  0.084  0.132  0.184 1.06   116
25      grIntT[3]  0.882 0.030  0.826  0.882  0.943 1.04    44
26      grIntT[4]  0.326 0.022  0.284  0.326  0.368 1.00   230
27      grIntT[5]  0.069 0.034 -0.005  0.069  0.135 1.40    59
28      grIntT[6]  0.169 0.033  0.109  0.168  0.236 1.06    78
29      grIntT[7]  1.649 0.048  1.562  1.647  1.751 1.01    54
30      grIntT[8]  0.723 0.054  0.616  0.724  0.825 1.01    42
31      grIntT[9] -0.167 0.082 -0.326 -0.163 -0.025 1.02    23
32     grIntT[10]  0.375 0.088  0.152  0.383  0.516 1.29    22
33     grIntT[11]  2.246 0.143  1.969  2.244  2.513 1.09    16
Code
  # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
  source('./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R')
  
  #summaryMod3_tt_growth_z0 <- MCMCSummaryRMNA(object = out_NN_OBmod3$mcmc, params = c("weight"), round = 3) %>%
  summaryMod3_tt_w <- MCMCSummaryRMNA(object = out_NN_OBmod3$mcmc, params = c("weight"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod3_tt_z <- MCMCSummaryRMNA(object = out_NN_OBmod3$mcmc, params = c("z"), round = 3) %>%
    rownames_to_column(var = "var") |> 
    mutate(
      ind0 = str_match(var, "([0-9]+), ([0-9]+)")[,2],
      t0 =   str_match(var, "([0-9]+), ([0-9]+)")[,3],
      ind = as.numeric(ind0),
      t = as.numeric(t0)
    ) |> 
    dplyr::select(-c(t0, ind0))
  
  summaryMod3_tt_z0 <- summaryMod3_tt_z |> 
    left_join(summaryMod3_tt_w, by = c("ind", "t"), suffix = c("_z", "_weight"))

Add other variables to summary values

Code
  ehLong <- 
    eh$eh |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      enc = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  firstLong <- eh$first |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(first = value)
  
  lastLong <- eh$last |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) |> 
    rename(last = value)
    
  cohortLong <- eh$cohort |> 
    as_tibble() |> 
    rownames_to_column("ind") |> 
    mutate(ind = as.numeric(ind)) 
  
  weightLong <- 
    eh$weight |>
    as.data.frame() |> 
    rownames_to_column("ind0") |> 
    pivot_longer(starts_with("ais")) |> 
    mutate(
      t = as.numeric(str_match(name, "([0-9]+)")[,1]),
      weightDATA = value,
      ind = as.numeric(ind0)
    ) |> 
    dplyr::select(-c(name, value, ind0))
  
  summaryMod3_tt_z <- summaryMod3_tt_z0 |> 
    left_join(ehLong) |> 
    mutate(
      # meanM1 = mean - 1,
      # pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
      enc8 = ifelse(enc == 1, 0.8, 0)
    ) |> 
    left_join(firstLong) |> 
    left_join(lastLong) |> 
    left_join(cohortLong) |> 
    left_join(weightLong) |> 
    mutate(
      resid = weightDATA - mean_weight
    ) 
    
  # summaryMod3_tt_z <- summaryMod3_tt_z0 |> 
  #   left_join(ehLong) |> 
  #   mutate(
  #     # meanM1 = mean - 1,
  #     # pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),
  #     enc8 = ifelse(enc == 1, 0.8, 0)
  #   ) |> 
  #   left_join(firstLong) |> 
  #   left_join(lastLong) |> 
  #   left_join(cohortLong) 
Code
  resids_mod3 <- summaryMod3_tt_z |> 
    group_by(ind) |> 
    summarize(
      meanResid = mean(abs(resid), na.rm = TRUE),
      n = n()
    ) |> 
    ungroup()

  summaryMod3_tt_zR <- 
    summaryMod3_tt_z |> 
    left_join(resids_mod3)
Code
ojs_define(numTags_OJS_mod3 = dim(eh$tags)[1])
ojs_define(summaryMod3_tt_zR_OJS = (summaryMod3_tt_zR))
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))

19.0.3.6 Plot predicted and observed masses for selected individuals

Select one or more individuals:

Code
viewof selectInd_mod3 = Inputs.select([...Array(numTags_OJS).keys()], {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})

summaryMod3_tt_zR_OJS_T = transpose(summaryMod3_tt_zR_OJS)
summaryMod3_tt_zR_OJS_T_selected = summaryMod3_tt_zR_OJS_T.filter((d) =>
  selectInd_mod3.includes(d.ind)
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code
Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Body mass (mg)" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    Plot.line(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_weight"
    }),
    <!-- Plot.dot(summaryMod3_tt_zR_OJS_T_selected, { -->
    <!--   x: "t", -->
    <!--   y: "mean_gr", -->
    <!--   fill: "green" -->
    <!-- }), -->
    Plot.dot(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "weightDATA",
      fill: "orange"
    }),
    Plot.ruleX(summaryMod3_tt_zR_OJS_T_selected, {
      x: "first",
      y: 3,
      "stroke": "green"
    }),
    Plot.ruleX(summaryMod3_tt_zR_OJS_T_selected, {
      x: "last",
      y: 3,
      "stroke": "red"
    }),
    Plot.text(summaryMod3_tt_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 60,
        frameAnchor: "top-left",
        text: (d) => "Cohort = " + d.cohort
      })
    ),
    Plot.text(summaryMod3_tt_zR_OJS_T_selected,
      Plot.selectFirst({
        x: 1,
        y: 56,
        frameAnchor: "top-left",
        text: (d) => "Residual = " + d.meanResid
      })
    )
  ],
  facet: {
    data: summaryMod3_tt_zR_OJS_T_selected,
    x: "ind"
  }
})

19.0.3.7 Plot survival

Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code
Plot.plot({
  width: width,
  height: 350,
  inset: 10,
  color: {
    scheme: "greys"
  },
  x: { label: "Age/season combination" },
  y: { label: "Probability of survival" },
  marks: [
    Plot.frame(),
    Plot.dot(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_z"
    }),
    Plot.dot(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "enc8",
      fill: "orange"
    }),
    Plot.line(summaryMod3_tt_zR_OJS_T_selected, {
      x: "t",
      y: "mean_z"
    }),
    Plot.ruleX(summaryMod3_tt_zR_OJS_T_selected, {
      x: "first",
      y: 1,
      "stroke": "green"
    }),
    Plot.ruleX(summaryMod3_tt_zR_OJS_T_selected, {
      x: "last",
      y: 1,
      "stroke": "red"
    })
    ,
    Plot.text(summaryMod3_tt_zR_OJS_T_selected, 
      Plot.selectFirst({
        x: 11,
        y: 1,
        text: d => d.cohort
      })
    )
  ],
  facet: {
    data: summaryMod3_tt_zR_OJS_T_selected,
    x: "ind"
  }
})

19.0.3.8 Output model summary data for Xiaowei

Code
  #write.csv(summary_tt, './models/cmrNN_OB/tt/dataOut/xiaowei/summary_tt.csv')
  
  #write.csv(summary_tt_z, './models/cmrNN_OB/tt/dataOut/xiaowei/summary_tt_z.csv')